#include "tbparameters.hpp"

void insertblocks(int & row, int & col,
		  cx_mat & H,
		  const cx_mat Md,
		  const cx_mat H1,
		  const cx_mat Mu
		  ){
  int dim=H1.n_rows;

  H.submat(row,col,row+dim-1,col+dim-1) = Md;
  col += dim;
  H.submat(row,col,row+dim-1,col+dim-1) = H1;
  col += dim;
  H.submat(row,col,row+dim-1,col+dim-1) = Mu;
  col -= dim;
  row += dim;
}


cx_mat transfermat(double E,
		   const cx_mat & H, 
		   const cx_mat & V1, 
		   const cx_mat & V2){
  int dim = H.n_rows;
  int dim2=2*dim;
  cx_mat T(dim2,dim2);
  cx_mat I(dim,dim);

  I.eye(dim,dim);
  T.zeros();
  
  T.submat(0,0,dim-1,dim-1)    = -solve(V1,H - E*I); 
  T.submat(0,dim,dim-1,dim2-1) = -solve(V1,V2.t());
  T.submat(dim,0,dim2-1,dim-1) = I;
  return T;
}
void superlatticeT(double E,
		   string name,string mat1, 
		   int nmat1,string mat2, 
		   int nmat2,rowvec kvec,
		   cx_mat & T1,cx_mat & T2, cx_mat & T12, cx_mat & T21){

  cx_mat  H1, H1c,  H1a, V1u, V1d;
  cx_mat  H2, H2c,  H2a, V2u, V2d;

  struct tb p1 = parameters(name,mat1);
  struct tb p2 = parameters(name,mat2);


  hamiltonian(p1,kvec,H1,H1c,H1a,V1u,V1d);
  hamiltonian(p2,kvec,H2,H2c,H2a,V2u,V2d);
  
  int dim = H1.n_rows;
  T1.zeros(2*dim,2*dim);
  T2.zeros(2*dim,2*dim);
  T12.zeros(2*dim,2*dim);
  T21.zeros(2*dim,2*dim);
  
  cx_mat X = 0.5*(V1u + V2u);
  cx_mat Y = X;
  cx_mat Ba = (H1a + H2a)*0.5;
  cx_mat Bc = (H1c + H2c)*0.5;

  int dim2 = 2*dim;
  
  T1  = transfermat(E,H1c,V1u,V1d) * transfermat(E,H1a,V1d,V1u);
  T2  = transfermat(E,H2c,V2u,V2d) * transfermat(E,H2a,V2d,V2u);
  
  T12 = transfermat(E,Bc,V2u,Y) * transfermat(E,Ba,Y,V1u);
  T21 = transfermat(E,Bc,V1u,X) * transfermat(E,Ba,X,V2u);
  
}

void superlatticeH(double E,
		   string name, 
		   string mat1, int nmat1,
		   string mat2, int nmat2,
		   rowvec kvec,cx_mat & H,
		   const bool velocitymat, const int cartesiancoord){

  cx_mat  H1, H1c,  H1a, V1u, V1d;
  cx_mat  H2, H2c,  H2a, V2u, V2d;


  struct tb p1 = parameters(name,mat1);
  struct tb p2 = parameters(name,mat2);


  if (not velocitymat){
    hamiltonian(p1,kvec,H1,H1c,H1a,V1u,V1d);
    hamiltonian(p2,kvec,H2,H2c,H2a,V2u,V2d);
  }else{
    velocitymatrix(p1,kvec,H1,H1c,H1a,V1u,V1d,cartesiancoord);
    velocitymatrix(p2,kvec,H2,H2c,H2a,V2u,V2d,cartesiancoord);
  }
  int dim1 = H1c.n_rows;
  int dim2 = H2c.n_rows;

  cx_mat X = 0.5*(V1d + V2d);
  cx_mat Y = X;
  cx_mat Ba = (H1a + H2a)*0.5;
  cx_mat Bc = (H1c + H2c)*0.5;
  
  int fdim = 2*(nmat1)*dim1 + 2*(nmat2)*dim2;
  int dim = dim1;
  H.zeros(fdim,fdim);

  cx_mat V1dt = V1d.t();
  cx_mat V1ut = V1u.t();
  cx_mat V2dt = V2d.t();
  cx_mat V2ut = V2u.t();
  cx_mat Ud   = (V1d+ V2d)*0.5;
  cx_mat Udt  = Ud.t();

  int dim4 = dim*4;


  H.submat(0,0,dim-1,dim-1)         = Bc;
  H.submat(0,dim,dim-1,2*dim-1)     = V1u;
  H.submat(0,fdim-dim,dim-1,fdim-1) = Udt;

  int row = dim;
  int col = 0;

  insertblocks(row,col,H,V1ut,H1a,V1d);
  
  for (int i=1;i<nmat1-1;i++){
    insertblocks(row,col,H,V1dt,H1c,V1u);
    insertblocks(row,col,H,V1ut,H1a,V1d);
  }

  insertblocks(row,col,H,V1dt,H1c,V1u);
  
  insertblocks(row,col,H,V1ut,Ba,Ud );
  insertblocks(row,col,H,Udt ,Bc,V2u);
  
  insertblocks(row,col,H,V2ut,H2a,V2d);

  for (int i=1;i<nmat2-1;i++){
    insertblocks(row,col,H,V2dt,H2c,V2u);
    insertblocks(row,col,H,V2ut,H2a,V2d);
  }

  insertblocks(row,col,H,V2dt,H2c,V2u);

  H.submat(fdim-dim,0,fdim-1,dim-1)                 = Ud;
  H.submat(fdim-dim,fdim-dim,fdim-1,fdim-1)         = Ba;
  H.submat(fdim-dim,fdim-2*dim,fdim-1,fdim-dim-1) = V2ut;

}


    /*
    H.submat(row,col,row+dim-1,col+dim-1) = V1dt;
    col += dim;
    H.submat(row,col,row+dim-1,col+dim-1) = H1c;
    col += dim;
    H.submat(row,col,row+dim-1,col+dim-1) = V1u;
    col -= dim;
    row += dim;
    H.submat(row,col,row+dim-1,col+dim-1) = V1ut;
    col += dim;
    H.submat(row,col,row+dim-1,col+dim-1) = H1a;
    col += dim;
    H.submat(row,col,row+dim-1,col+dim-1) = V1d;

    col -= dim;
    row += dim;
    */
